function A407()
format long;

% 基于密度的空间聚类（density-based spatial clustering of applications with noise，DBSCAN）
% DBSCAN算法将样本观测点视为聚类变量空间中的点，以样本观测点O邻域内的邻居个数作为O所在区域的密度测度。

% === 读取保存的 .mat 数据 ===
load('gprdata.mat');  % 包含 data, coords, rh_nsamp, scansnum

% === 聚类数据裁剪与Z轴反转 ===
X = coords;
xMin = 1; xMax = size(data,2);
yTarget = 0; zMin = 1; zMax = 512;
indices = X(:,1) >= xMin & X(:,1) <= xMax & ...
          X(:,2) == yTarget & ...
          X(:,3) >= zMin & X(:,3) <= zMax;
X_cropped = X(indices, :);
X_cropped(:,3) = max(X_cropped(:,3)) - X_cropped(:,3);  % 反转Z

subplot(2, 2, 1)
scatter(X_cropped(:,1), X_cropped(:,3), '.');
xlabel('X坐标'); ylabel('反转后的Z坐标');
title('原始数据 (Z轴反转展示)');
xlim([xMin xMax]);
ylim([0 max(X_cropped(:,3))]);
grid on;

% === 特征构建与标准化 ===
weight_x = 0.01;
weight_z = 1;
X_features = [X_cropped(:,1)*weight_x, X_cropped(:,3)*weight_z];
X_features = rescale(X_features);

% === 自动选择 epsilon（最大曲率法）===
minpts = 5000;
[~, D] = knnsearch(X_features, X_features, 'K', minpts+1, 'Distance', 'cityblock');
k_distances = D(:, end);
sorted_k_dist = sort(k_distances);

nPoints = length(sorted_k_dist);
allPts = [1:nPoints; sorted_k_dist']';
firstPoint = allPts(1,:);
lastPoint = allPts(end,:);
lineVec = lastPoint - firstPoint;
lineVecNorm = lineVec / norm(lineVec);
vecFromFirst = allPts - firstPoint;
scalarProduct = vecFromFirst * lineVecNorm';
vecProjected = scalarProduct * lineVecNorm;
vecToLine = vecFromFirst - vecProjected;
distToLine = sqrt(sum(vecToLine.^2, 2));
[~, idx] = max(distToLine);
epsilon = sorted_k_dist(idx);

disp(['选择的 minpts: ', num2str(minpts)]);
disp(['选择的 epsilon: ', num2str(epsilon)]);

subplot(2, 2, 2)
plot(sorted_k_dist);
hold on; plot(idx, epsilon, 'ro');
title('k-distance图 (cityblock距离)');
xlabel('点排序'); ylabel([num2str(minpts) '个最近邻距离']);
grid on;

% === DBSCAN 聚类并可视化 ===
labels1 = dbscan(X_features, epsilon, minpts, 'Distance', 'cityblock');
subplot(2, 2, 3)
if any(labels1 == -1)
    X_plot = X_features(labels1 ~= -1, :);
    labels_plot = labels1(labels1 ~= -1);
    gscatter(X_plot(:,1), X_plot(:,2), labels_plot);
    title(['DBSCAN聚类 (epsilon=', num2str(epsilon), ', minpts=', num2str(minpts), ')']);
else
    gscatter(X_features(:,1), X_features(:,2), labels1);
    title(['DBSCAN聚类 (epsilon=', num2str(epsilon), ', minpts=', num2str(minpts), ')']);
end
xlabel('标准化X坐标'); ylabel('标准化Z坐标'); grid on;

% === 绘制灰度图 ===
subplot(2, 2, 4)
sampling_interval = 25 / rh_nsamp;
Y = sampling_interval:sampling_interval:25;
xmax = (scansnum-1) * 0.01;
X_axis = 0:0.01:xmax;

imagesc(X_axis, Y, data);
set(gca, 'YDir', 'reverse');
axis([0 xmax min(Y) max(Y)]);
colormap(gray);
title('B-scan灰度图');
xlabel('里程 x/m'); ylabel('时间轴 t/ns');

end
